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CHARACTERIZATION  OF  NONLINEAR  SYSTEMS  WITH  MEMORY 
BY  MEANS  OF  VOLTERRA  EXPANSIONS  WITH  FREQUENCY  PARTITIONING: 
APPLICATION  TO  A  CICADA  MATING  CALL 


1.  INTRODUCTION 


The  Wiener  characterization  of  time-invariant  nonlinear  systems  with  memory  utilizes  an 
expansion  involving  kernels  of  various  orders,  typically  limited  to  third  order  or  less  in  practice. 
This  limitation  to  third  order  is  due  to  the  “curse  of  dimensionality,”  that  is,  the  fact  that  the 
required  number  of  kernel  coefficients  increases  exponentially  with  the  order  of  fit  adopted  in  the 
expansion.  Also,  great  care  must  be  taken  in  keeping  the  condition  number  of  the  very  large 
relevant  data  matrix  within  reasonable  bounds,  so  that  its  pseudo-inverse  and  Wiener  expansion 
are  accurate.  This  must  be  achieved  even  when  the  Gaussian  excitation  has  a  colored  spectrum 
and.  therefore,  requires  a  special  type  of  expansion  of  the  kernel  terms  themselves.  The  end 
result  of  the  characterization  is  a  quantitative  statement  and  breakdown  of  the  exact  amounts  of 
each  type  of  nonlinearity  at  the  nonlinear  system  output.  This  information  has  obvious 
applications,  including  the  underwater  acoustic  channel,  by  indicating  the  strength  of  the  system 
nonlinearities  and  whether  they  can  be  ignored  or  whether  they  should  be  exploited  for  improved 
detection  and/or  classification  purposes. 

However,  there  are  numerous  instances  where  the  excitation  .v(/)  of  the  nonlinear  system 
under  investigation  is  not  under  one’s  control,  but  is  instead  dictated  by  the  particular  physical 
situation  at  hand.  Such  a  situation  arises  in  the  case  of  a  physical  cicada  mating  call,  where  the 
excitation  is  the  mechanical  movement  of  the  ereature’s  tymbal(s),  and  the  nonlinear  system 
response  z(t)  is  the  audio  pressure  waveform  received  at  a  recording  device  some  distance  away. 
In  such  a  ease,  appeal  cannot  be  made  to  the  Wiener  expansion  because  the  higher  order 
moments  of  excitation  x(t)  are  then  unknown.  Whereas  all  the  odd-order  moments  of  a 
stationary,  zero-mean,  Gaussian  proeess  are  zero,  and  all  the  even-order  moments  arc  expressible 
in  terms  of  just  the  correlation  function  Rx(t)  of  the  Gaussian  process  x(t),  none  of  the  higher 
order  moments  of  a  non-Gaussian  process  are  generally  known,  and  arc  extremely  difficult  to 
estimate  due  to  the  need  for  an  extreme  amount  of  data  about  the  excitation  process  x(t).  For  this 
reason,  in  the  current  cicada  application,  it  is  necessary  to  abandon  the  Wiener  procedure  for 
third  order  and  above  and,  instead,  revert  to  the  Voltcrra  expansion,  which  makes  no 
assumptions  about  the  statistics  of  the  excitation  x(t). 

In  an  attempt  to  alleviate  the  curse  of  dimensionality,  a  frequency  partitioning  scheme  was 
proposed  in  an  earlier  technical  report.*  However,  it  has  been  found  that  the  non-overlapping 
square  regions  there,  in  the  two-dimensional  frequency  space  fufi,  interact  with  each  other, 
cross-contaminating  each  estimate.  Also,  off-diagonal  frequency  squares  bring  up  the  question 


*Albert  H.  Nuttall,  “Characterization  of  Nonlinear  Systems  with  Memory,  by  Means  of  Volterra  Expansions,” 
Adaptive  Methods  Inc.,  Middletown,  RI,  30  December  2009. 


as  to  what  frequency  band  to  fit,  the  lower  or  upper  or  both.  These  problems  have  led  to  a  very 
different  partitioning  of  the  two-dimensional  frequency  space,  to  be  described  here,  that  solves 
both  problems  and  yields  kernel  and  waveform  estimates  of  both  first  and  second  order  that  are 
excellent  for  the  control  example  used  to  test  this  new  technique.  Numerical  results  are 
presented  to  demonstrate  the  efficiency  of  this  latest  approach.  An  application  to  a  cicada  mating 
call  is  made  that  illustrates  a  very  interesting  second-order  behavior. 
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2.  BASIC  VOLTERRA  MODEL 


Nonlinear  system  excitation  x(t)  is  sampled  at  frequency/*  Hz,  resulting  in  time-sampling 
increment  A  =  1  lfs  seconds  and  sampled  sequence  {.v(/7  A)}.  For  simplicity  of  notation,  the  A 
symbol  will  frequently  be  suppressed  and  the  excitation  sequence  will  be  denoted  simply  by 
{.v(/7)}.  At  other  times,  A  will  be  kept  in  order  to  stress  the  time  dependence.  Also,  sequence 
{.v(/7)}  will  frequently  be  referred  to  as  a  waveform. 

Consider  a  time-invariant  nonlinear  system  with  actual  sampled  input  sequence  {x(/?)}  and 
actual  sampled  output  sequence  {z(/7)j ,  both  of  w  hich  are  sampled  at  the  same  rate  fs  and 
recorded  simultaneously.  The  causal  time-invariant  Volterra  mode /  sampled  output  sequence 
{ y(n)}  is  then  given,  to  third  order,  by 

y(/7)  =  hQ  +  £  /?,  (£, )  x( n  -  k{ )  +  £  £  h2  (k , , k2 )  x(n  -  A, )  x{n  -  k2 ) 

jt,=0  *,=0  k2=() 

+  Z  X  x(n-k2)  x(n-k ,)  (1) 

Aj=0  A:=0  A =0 

=  To  +  y\  (n) + y2  (n)  +  y ,  (/?), 

where  h(t ,  ht ,  h2 ,  /7,  are  the  zeroth-order  through  third-order  (time-invariant)  time-domain  kernels 
of  the  Volterra  expansion.  It  is  assumed  that  the  Volterra  kernels  /?, ,  h2 ,  hy  are  represented  with 

the  same  time-sampling  increment  A  as  used  for  the  nonlinear  system  input  and  output 
waveforms  x(n)  and  z(n).  (The  {  }  symbol  around  a  sequence  will  be  dropped  from  here  on.) 

It  is  also  assumed  that  the  same  “memory  length”  K  in  equation  (1 )  is  appropriate  for  all  three 
orders  of  these  kernels.  This  may  not  be  the  case  in  practice;  if  not,  different  sizes  K{,  A',,  A',  of 

the  summations  in  equation  (1)  may  be  required  for  a  decent  approximation  of  nonlinear  system 
output  :(n)  by  the  model  output  y(/7). 

Although  the  input  data  x(n)  appear  in  a  nonlinear  fashion  in  component  model  outputs 
y2(n)  and  y}(n)  in  equation  (1),  sequence  x(/?)  is  known ,  and  the  required  second-order  and 
third-order  input-data  products  in  equation  (1 )  can  be  easily  calculated  and  stored.  The 
unknowns  in  the  Volterra  expansion  (1)  are  the  four  kernels  h() ,  /?, ,  h2 ,  h2 ,  which  appear  linearly 

in  the  model  output  v(/7).  This  observ  ation  strongly  suggests  the  use  of  a  least  squares  approach 
in  attempting  to  fit  model  output  y(n)  to  the  actual  measured  nonlinear  system  output  z(n), 
through  optimum  choice  of  these  Volterra  kernels;  see  figure  1.  Namely,  the  equations 
determining  the  best  kernels  ha ,  /?, ,  h2 ,  h2  will  be  the  solutions  of  simultaneous  linear  equations 
under  the  least  squares  philosophy.  This  is  a  very  important  consideration. 
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Figure  1.  Measurement  and  Fitting  Procedure 


One  major  problem  with  this  Volterra  approach  is  the  ill-conditioning  of  the  very  large  data 
matrix  that  must  be  pseudo-inverted  in  the  least  squares  processing  approach.  For  a  Gaussian 
excitation  x(/7),  this  ill-conditioning  problem  could  be  partially  overcome  by  switching  to  a 
Wiener  expansion  with  uncorrelated  components  (on  an  ensemble  average  basis);  but  that  route 
is  not  available  for  third  order  and  above  in  the  current  case  of  a  non-Gaussian  excitation  jc(/?). 
However,  a  second-order  Wiener  modification  is  still  applicable  and  is  described  below. 

The  additional  problem  of  very  large  data  matrices  can  be  overcome  by  time-partitioning  of 
the  excitation  x(n),  accompanied  by  temporary  storage  on  a  hard  drive.  When  particular  sub¬ 
matrices  are  needed  for  computation  in  random  access  memory  (RAM),  these  matrices  can  be 
recalled  as  needed  from  the  hard  drive.  Of  course,  the  tradeoff  here  is  that  the  execution  time 
can  increase  considerably. 

However,  the  really  major  problem  associated  with  both  the  Volterra  and  the  Wiener 
expansions  is  the  curse  of  dimensionality  (COD),  namely,  the  extreme  number  of  coefficients 
(kernel  values)  required  in  equation  ( 1 ).  At  first  order,  the  number  of  coefficients  that  must  be 
determined  is  A/,  =  AT;  at  second  order,  the  number  of  coefficients  is  approximately 

M2  =  K2 /l ;  and  at  third  order,  it  is  approximately  A/,  =  K' /b .  The  helpful  denominator 
factors  of  2  and  6  come  about  through  the  judicious  imposition  of  symmetry  on  the  unknown 
kernels  /?,  and  /?3,  without  loss  of  generality. 

The  programming  required  for  these  tasks  requires  careful  consideration  in  order  to  avoid 
ill-conditioning,  storage  limitations,  and  excessive  execution  time.  Once  available,  the  resulting 
program(s)  can  be  utilized  to  determine  the  amounts  and  orders  of  nonlinearity  in  the  cicada 
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mating  call.  If  there  is  significant  energy  in  the  higher  order  nonlinear  terms,  this  will  indicate 
quantitatively  just  how  important  it  is  to  determine  exactly  how  the  cicada  can  emit  such  strong 
signals  from  its  small  body.  Practical  applications  of  this  information  are  of  obvious  importance 
in  many  fields,  including  underwater  acoustics. 

For  an  expansion  limited  to  second  order,  the  Wiener  modification  takes  the  form 
>-(/?)  =  u’0  +  £  vv, (£, )  x(n - k, )  +  £  vc, (k, ,k2 )  [v( n - kx )  x(n -k2)~  Rx {k{  - k2 )] 

k  |  =0  k |  =0  k j  =0 


=  “  Z  Z  W2  (k\  ’  k2  )  RM\  -  k2  ) 

At=0  A-:=0 
A'  l  AT  - 1 

+  Z  Z  x(n-k2) 

Al=  0  k:~0 

=  y0 +>',(«)+ v2(//), 

w  here  w0 ,  vv, ,  vv,  are  the  Wiener  kernels.  The  subtraction  of  the  excitation  correlation  function 
(k{  -  k , )  in  the  first  line  of  equation  (2)  makes  the  second-order  (double-summation)  term  in 
the  top  line  uncorrelated  with  the  zeroth-order  term  vrn ,  regardless  of  the  statistics  of  excitation 
x(n).  If  x(n)  has  zero  mean,  then  the  zeroth-order  term  and  the  first-order  term  in  the  first  line 
of  equation  (2)  are  also  uncorrelated  w  ith  each  other.  But  that  still  leaves  the  cross-correlation 
between  the  first-order  term  and  the  second-order  term  unknown  in  the  case  of  a  non-Gaussian 
excitation  x(n).  Nevertheless,  this  Wiener  modification  in  equation  (2)  is  well  worth 
considering  and  trying,  when  the  expansion  is  limited  to  second  order.  It  can  lead  to 
significantly  smaller  condition  numbers  for  the  large  data  matrix  involved  in  the  least  squares 
procedure.  Only  the  correlation  function  Rx(k)  of  the  excitation  jc(/?)  needs  to  be  estimated 
from  the  available  data;  this  is  a  very  simple  and  quickly  done  task. 

The  interrelationships  between  the  Volterra  and  Wiener  kernels  at  second  order  are  seen, 
from  equations  ( 1 )  and  (2),  to  be 

k  \  K  \ 

K  =  w’o  -  Z  Z  w2  (A| ,k2)  Rx{k ,  -k2),  /7,  =  vv,,  /?,  =  vv, .  (3) 

*1=0  *,=0 

Thus,  if  the  first  line  of  Wiener  model  (2)  is  used  to  least  squares  fit  the  measured  nonlinear 
output  -(/?),  the  corresponding  Volterra  kernels  can  then  be  determined  from  equation  (3),  and 
the  corresponding  Volterra  waveforms  can  be  determined  from  the  last  two  lines  of  equation  (2). 


*1=0 
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For  completeness,  the  Wiener  modification  at  third  order  for  a  Gaussian  excitation  would  be 

Af-1  K  1  K- 1 

y(n)  =  w '0  +  £  W|(*,)jr(i*-A,)  +  Z  ^  w2(*„*2)[jf(#i-^l)x(ii-*2)-J?I(*1  -*2)] 

*i=0  A'j=0  f<2  —  0 

(4) 

+  Z  Z  Z  W«-A,)  x(n—k2)  x(n-ki)-3  Rx{k2  -  k2)  x(n  -  A,)], 

A,  =0  k2=  0  *3=0 

using  the  assumed  symmetry  of  Wiener  kernel  wy(kl9k29ky)9  without  loss  of  generality.  The 
corresponding  Volterra  kernels  are  then  given  by 

K  l  A'  -  1 

/?0  —  W0  —  ^2  (^1^2)  *,(*,  ”^2)’ 

A,  =0  *,=0 

Ai (A, )  =  w, (A, )  —  3  Y,  Z  w’i(ApA2»Aj)  /?,(*, -A,),  (5) 

A2=0  *,=0 

h2  =  u’: ,  h2  =  vc, , 

and  would  be  computed  after  the  least  squares  fit  of  Wiener  model  (4)  to  measured  nonlinear 
system  output  data  z(n)  was  completed.  Volterra  waveforms  >’0  ■>  >’1  V2  (")» >’3  (/7)  are  still 
given  by  equation  (1)  in  terms  of  the  Volterra  kernels. 
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3.  FIRST-ORDER  CONSIDERATIONS 


3.1  TIME-BANDWIDTH  PRODUCT  PROBLEM 

Corresponding  to  the  first-order  time-domain  description  /?,(A)  =  /?,(£, A)  in  equation  (1),  it 
is  useful  to  consider  the  equivalent  frequency-domain  descriptor,  namely,  the  first-order 
frequency-domain  kernel  (discrete-time  Fourier  transform) 

//,(/,)  =  !  /?,(A:1A)exp(-/2^AlA/1).  (6) 

*,=o 

This  relation  can  be  inverted  by  means  of  the  inverse  Fourier  transform  to  obtain 

F 

/?,(A,A)  =  A  Jr//',  //,(/,)  exp(/2^A,A/1),  for  A,  =  0:A-1,  (7) 

F 

where  the  assumed  band-limited  character  of  the  frequency-domain  kernel  //,  (/, )  has  been 
explicitly  indicated  as  (~F,F).  That  is,  F<  l/(2A)  =  /s/2,  the  Nyquist  frequency. 

Suppose  that  first-order  time-domain  kernel  /?,  (A,  A)  in  equations  ( 1 )  and  (6)  has  a  time 
extent  of  T  seconds, 

T  =  K  A,  (8) 

and  that  the  first-order  frequency-domain  kernel  //,  (/, )  in  equation  (6)  has  a  bandwidth 
(frequency  extent)  of  F  Hz,  assumed  lowpass  about  zero  frequency,  as  indicated  in  equation  (7). 
Then,  to  get  a  complete  characterization  of  the  continuous  function  /7,(r),  it  is  necessary  to 
sample  ht  (r)  with  time  increment  1/(2 F)  seconds  or  finer.  Since  the  time  duration  of  /?,(r)  is  T 
seconds,  a  total  of  at  least  2TF  samples  must  be  employed  in  order  not  to  miss  any  significant 
information  about  A,(r).  This  implies  that  the  size  K  in  equation  (1)  must  be  at  least  2 TF. 

Although  this  value  K  =  2 TF  is  usually  manageable  at  first  order,  as  far  as  least  squares  is 
concerned,  the  numbers  of  required  coefficients  for  specification  at  second  and  third  orders, 
given  by 


M2=(2TF)2/2  and  A/,  =(27F)3/6,  (9) 

respectively,  can  quickly  get  out  of  hand.  For  example,  if  2 TF  =  100,  then 

A/,  =  A  =  100,  M2  =5,000,  and  A/,  =167,000.  (10) 
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In  the  normal  equations  that  arise  in  least  squares,  the  size  of  the  data  product  matrix  that  must 
be  inverted  is  M  x  A/.  The  M2  x  A/,  case  can  often  be  solved  with  current-day  computer  RAM, 
but  the  A/,  x  A/,  matrix  will  often  not  fit  into  RAM.  If  a  simultaneous  fit  of  all  the  components 

in  equation  (1)  to  measured  nonlinear  system  output  z(n)  were  of  interest,  it  can  happen  that 
nothing  at  all  could  be  achieved,  because  of  excessive  storage  requirements.  Some  method  of 
addressing  partitions  of  the  various  kernels  is  needed  if  meaningful  useful  estimates  of  the 
various  kernels  are  to  be  obtained  at  higher  orders. 


3.2  ALLEVIATION  OF  THE  TF  PRODUCT  AT  FIRST  ORDER 

To  see  how  a  large  TF  product  can  be  alleviated,  recall  equation  (7).  For  a  real  first-order 
kernel  A,  (A, A),  develop  it  as 


A, (A, A)  =  2 A  Re  //,(/;)  exp(/ 2  ;r  A,  A  /, ) 

0 

IV  21V 

=  2A  Re  J z//,  )  exp(/2^A-,A /,)  +  2A  Re  //,(/,) exp(i2/r A/,) +  •••  (1 1) 

o  w 

=  A10(A,A)  +  A, ,  (A, A)  +  h[2  (A, A)  +  •  •  •  +  hu  (A, A),  J  =  F/W- 1 . 


Subscript  1/  denotes  the  j-th  band  of  the  first-order  kernel.  This  breakdown  of  kernel  A,  (A,  A)  is 
not  an  approximation;  it  is  exact.  The  frequency  content  of  component  kernel  A10(A,A),  namely, 
its  frequency-domain  kernel  Hw (/,),  is  limited  to  positive  frequency  band  (0,  IT),  while 
//,-(/,)  for  component  kernel  A,  (A,  A)  is  limited  to  positive  frequency  band  (jW  ,(j  +\)W)  for 
j  -  0 :  J .  This  reduces  the  time-bandwidth  product  for  each  component  kernel  by  a  factor  of 
F/W .  On  the  other  hand,  there  are  now  J  + 1  =  F/IV  narrowband  kernels  to  be  determined 
instead  of  just  one  broadband  kernel.  If  one  can  evaluate  (estimate)  all  the  individual  kernels 
A,  (A, A),  j  =  0  :  J,  they  can  be  added  together  to  yield  the  total  first-order  time-domain  kernel 

A,  (A, A),  and  thereby  determine  the  entire  first-order  model  output  y^n)  in  equation  (1). 

The  first-order  model  output  sequence  >’,(»)  in  equation  (1)  can  now  be  expressed  as  the 
sum  ofJ+  1  narrowband  components: 


y,  (w)  =  v10  (n)  +  ,(«)  +  •••  +  yu  (n) , 


(12) 


where 


Ti7(")  =  X  x(n~k\)  =A  J  d.f\  H\j(f\)  X  (/, )  exp(/  2  n  j\  n  A) , 

*,=0 


(13) 
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and  where 


X(f)  =  Yj  exp(-/2^/wA)  x(nA) 


(14) 


n 


is  the  Fourier  transform  of  the  discrete-time  excitation  x(nA).  Regardless  of  the  spectral  content 
X(f)  of  time-excitation  x{nA),  a  band-limited,  first-order,  frequency-domain  kernel  //,,(/) 
w  ill  limit  the  spectral  content  of  model  output  component  y, ,  (n)  to  the  frequency  band 
(jW ,(  /  +  1)IF).  If  this  particular  model  component  y,  An)  is  to  be  separately  fitted  to  the 

measured  system  output  z(n),  it  makes  sense  to  first  filter  the  nonlinear  system  output  z(n)  to  this 
same  frequency  band  under  consideration.  Call  this  filtered  version  z  fin)  for  the  y'-th  band, 

y  =  0:  J. 

If  this  pre-filtering  of  z(n)  is  not  conducted,  the  least  squares  fit  of  bandpass  sequence  yl;.(«) 

to  the  original  nonlinear  system  output  sequence  z(n)  will  result  in  the  least  squares  procedure 
attempting  to  match  out-of-band  frequency  components,  which  will  be  considered  as  noise  or 
outliers  to  be  fitted  as  well  as  possible.  Therefore,  it  is  mandatory  to  remove  these  contaminating 
out-of-band  components  in  z(n)  before  attempting  a  least  squares  fit  of  y,  ( n )  to  the  available 

system  output  data.  This  eliminates  out-of-band  frequency  components  in  z(n)  that  the  model 
components  can't  possibly  fit  because  of  their  intentionally  limited  frequency  content.  Even  if  this 
partitioning  of  the  frequency  scale  were  not  adopted  at  first  order,  nonlinearity  output  process  z(n) 
should  still  be  filtered  to  the  model’s  band  of  interest  before  doing  any  least  squares  fitting, 
because  a  model  can’t  match  frequency  components  outside  its  own  band. 

The  breakdown  of  the  first-order  time-domain  kernel  in  equation  (1 1)  is  usually  not 
necessary  in  practice  because  the  TF  product  is  typically  not  that  large,  and  the  corresponding 
least  squares  data  matrix  can  be  easily  fitted  into  RAM.  The  real  benefit  of  reducing  the  TF 
product  is  realized  much  more  strongly  when  attempting  second-  and  third-order  least  squares 


fits. 


3.3  SINE  AND  COSINE  EXPANSION 

For  the  first-order  kernel  (waveform)  /?l0(r)  with  duration  T and  positive-frequency  extent 

(0 fiV),  an  attractive  expansion  (including  negative  frequencies)  is  afforded  by  the  standard 
trigonometric  form 


(15) 


where  M/T  =  IV,  the  highest  frequency  in  this  particular  low-frequency  band,  and  the  {a(m)} 
and  {h(m)\  coefficient  sequences  are  real.  The  number  of  coefficients  to  be  determined  is  about 
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2M- 2TW,  not  2TF.  Thus,  the  frequency-partitioning  procedure  achieves  a  worthwhile 
reduction  in  the  sizes  of  the  individual  fitting  procedures.  There  is  no  b( 0)  term  in  this  usual 
harmonic  expansion  in  equation  ( 1 5);  this  is  a  trivial  point  at  first-order,  but  will  become 
significant  at  second-order.  Similar  harmonic  expansions  to  equation  (15)  hold  for  the  frequency 
bands  (W,2W),  (2W -W ,F). 

The  first-order  time-domain  kernel  hi(kl A)  is  forced  to  be  real.  Therefore,  the  first-order 
frequency-domain  kernel  //,  (/J )  defined  in  equation  (6)  satisfies  a  conjugate  symmetry  property: 

//,(-/■)  =  (/,)*•  (16) 

If  one  specifies  //,  (/, )  for  positive  frequencies  /, ,  its  values  for  negative  frequencies  are 
automatically  fixed.  This  positive  frequency  range  is  called  the  “fundamental  region”  at  first 
order. 

A  single  positive-frequency  elemental  component  of  hl(kiA)  in  this  fundamental  region  is 


m 


rn . 


[a(m)  -  i  b{m)]  exp  /  2  n  —  A  ,  a(m),  b(m)  real;  fm= — >  0. 


(17) 


Upon  addition  of  the  conjugate  contribution  from  the  corresponding  negative  frequency  -m/T , 
the  fundamental  real  basis  function  for  this  single  frequency  fm  is,  after  scaling  by  1/2, 


m 


a(/«)  cos  \2n — A:,  A  +b{m)  sin 

\  T  ) 


m  ^ 
271— -k,  A 


J 


(18) 


This  is,  of  course,  the  same  expression  as  already  given  in  equation  (15);  however,  this  latter 
development  serves  to  introduce  the  “fundamental  region”  concept,  which  will  take  on  a  much 
greater  significance  at  second  and  higher  orders.  The  basis  functions  are  simply  cosines  and 
sines  at  first  order. 


3.4  KEY  OBSERVATION  BY  MEANS  OF  THE  FREQUENCY  DOMAIN 

The  first-order  component  output  time  waveform  _y.  (nA)  in  equations  (12)  and  (13)  was 
expressed  in  terms  of  frequency-domain  quantities  as 

>’iy  (,?A)  =  A  {  df  exp(/  2n  fn  A)  (19) 

where  the  excitation  spectrum  X(  f)  is  generally  broadband.  The  only  place  where  time 
variable  nA  appears  on  the  right-hand  side  of  this  expression  is  inside  the  exp(  )  with  the 
multiplier  f  If  waveform  v, ,  (nA)  on  the  left-hand  side  of  equation  ( 1 9)  is  to  contain  frequency 
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components  only  in  frequency  band  (jW,(j  + 1) IT),  then  frequency-domain  kernel  HXj(f )  must 
be  nonzero  only  in  this  same  band.  That  is,  one  must  have  //,,(/)  nonzero  only  for 

jW<f<U  +  \)W  (20) 

and  the  corresponding  negative  frequencies.  This  appears  to  be  a  trivial  point  here  at  first  order; 
however,  this  observation  will  be  of  crucial  importance  at  second  and  higher  orders. 


3.5  FIRST-ORDER  WAVEFORM  EXPANSION 

When  the  trigonometric  expansion  for  /?in(A',A)  in  equation  (15)  is  substituted  into  the  first 
equality  on  the  left-hand  side  of  equation  (13),  and  the  summations  are  interchanged,  the  result  is 

M  M 

y10(??A)  =  Yj  a(m)  xAn A,w)  +  £  b(m)  xs(nA ,»/),  (21) 

m-0  m  1 


where 


xc(nA,m) 

x^(nA,m) 


i\  —  1 

I 

K  I 

I 


cos 


27r  —  kA 


sin 


T 

m 

T 


2n—kA 


x(nA  -  kA) , 
x(nA  -  k A). 


(22) 


These  latter  summations  are  convolutions  of  the  trigonometric  basis  functions  w  ith  the  excitation 
.y(hA).  These  waveform  sequences  { jc(. (/?A, m) }  and  {av(/iA,w)}  can  be  precomputed  for  any  m 
values  of  interest  and  stored.  According  to  the  summation  limit  in  equation  (2 1 ),  these 
waveforms  are  approximately  band  limited  to  (0,1^  Hz.  as  desired,  since  M/T  =  W .  Expansions 
similar  to  equation  (2 1 )  can  be  developed  for  the  other  frequency  bands  (IT,  21V),...,(F -  W,  F). 

The  expansion  of  first-order  component  waveform  yl0(nA)  in  equation  (21)  is  linear  in  the 
unknown  real  coefficients  a(in)  and  b(m).  These  latter  coefficients  must  now  be  chosen  so  that 
v„,(»A)  approximates  z()(/?A)  as  closely  as  possible  in  a  least  squares  sense;  sec  the  discussion 

following  equation  (13).  The  result  of  the  least  squares  minimization  will  be  a  set  of  2M-  1 
simultaneous  linear  equations  for  the  unknown  real  coefficients.  (There  will  be  2M  simultaneous 
linear  equations  for  the  real  coefficients  in  all  the  other  bands,  j  =  1 :  J.) 

Once  the  coefficients  a(m)  and  b(m)  have  been  determined  from  the  least  squares  fit  to 
rn(nA),  equation  (15)  is  used  to  find  the  estimated  kernel  /710(A,A).  Also,  equations  (21)  and 
(22)  then  yield  the  best-fitting  waveform  j’l0(nA)  for  the  (0,lf)  frequency  band.  This  fitting 
procedure  must  be  repeated  for  all  the  other  IT-wide  bands  up  to  frequency  F. 


After  accomplishing  the  fits  for  all  the  frequency  bands  up  to  frequency  F,  producing  fitted 
first-order  waveforms  >’10 («A),  yn(nA),...,yu(nA),  these  waveforms  are  added  together  to  get 

the  total  first-order  fit  to  the  total  nonlinear  broadband  measured  response  z(n A).  Then, 

a  total  error  waveform  at  first  order  can  be  computed  according  to 

e,(>7A)  =  >’,(>7A)-z(«A).  (23) 


3.6  ELIMINATION  OF  ERRORS  AT  FREQUENCY  JOINTS 

It  has  been  found  that  the  spectrum  of  error  e]  (nA)  in  equation  (23)  has  peaks  in  the 
neighborhoods  of  frequencies  W,  2W,31V,  ....  These  latter  frequencies  are  the  locations  of  the 
“joints”  of  the  frequency  partitioning  in  equation  (11).  The  reason  for  the  errors  at  these 
frequency  joints  is  the  inability  to  do  perfect  bandpass  filtering  on  output  process  z(nA).  Every 
physical  filter  has  to  have  a  transition  region  in  frequency,  from  its  full  response  to  its  maximum 
attenuation  level.  Also,  the  finite-duration  (gated)  sines  and  cosines  employed  in  expansion  (15) 
always  have  side  lobes  in  the  frequency  domain  and  will  therefore  have  some  nonzero  response 
to  frequency  components  outside  their  nominal  “ideal”  band. 

A  method  for  overcoming  this  problem  is  best  explained  by  taking  a  numerical  example. 

Let  bandwidth  W  -  6  kHz.  First,  perform  a  least  squares  fit  over  (0:6)  kHz.  Retain  the 
coefficients  0(777)  and  6(777)  only  corresponding  to  the  band  (0:5)  kHz.  Discard  the  edge 
coefficients  corresponding  to  (5:6)  kHz,  because  they  are  expected  to  be  in  error. 

Next,  perform  a  separate  fit  over  the  (4:10)-kHz  band.  Retain  only  the  coefficients 
pertaining  to  the  interior  (5:9)-kHz  band.  Discard  both  edge  coefficients  corresponding  to  bands 
(4:5)  kHz  and  (9: 10)  kHz.  Next,  fit  over  (8:14)  kHz.  Retain  only  the  coefficients  for  (9:13)  kHz. 
The  general  procedure  is  now  obvious. 

In  this  manner,  although  6  kHz  is  fit  at  a  time,  only  the  interior  4-kHz  coefficients  are 
retained,  because  these  are  the  coefficients  expected  to  be  correct.  The  overlaps  in  frequency,  for 
each  successive  fit,  are  necessary  in  order  to  be  able  to  discard  edge  coefficients  that  are  not 
considered  trustworthy. 


3.7  FILTERING  THE  NONLINEAR  RESPONSE 

By  design,  component  waveform  jih(/7A)  has  frequency  content  only  in  the  frequency  band 
(O.JT);  of  course,  this  is  only  approximately  true,  as  noted  above.  But  nonlinear  system  output 
z(7?A)  is  very  broadband.  As  discussed  in  the  sequel  to  equation  (13),  model  output  y10(77A) 

cannot  possibly  match  the  frequency  components  of  z(nA)  that  are  outside  the  band  (0,JT), 
whether  they  are  signal  or  noise  components.  On  the  other  hand,  least  squares  is  a  greedy 
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procedure  and  will  inherently  attempt  to  fit  out-of-band  components  to  some  extent.  Therefore, 
filtering  of  z(nA)  is  necessary  before  attempting  any  least  squares  fits. 

In  doing  the  bandpass  filtering  of  nonlinearity  response  z(nA)  to  the  Ik-wide  band  of  the 
current  (j-th)  fit,  no  distortion  of  the  in-band  components  of  Zj(nA)  can  be  allowed  to  occur. 
Therefore,  it  is  necessary  to  utilize  a  bandpass  filter  with  a  linear  phase  shift  over  the  central  4 
kHz  of  each  sub-band  of  width  6  kHz.  A  symmetric  finite-impulse-response  (FIR)  filter  fits  this 
requirement  perfectly.  In  fact,  the  phase  response  of  such  a  filter  is  linear  for  all  frequencies. 
Also,  the  amplitude  response  of  each  filter  must  be  extremely  flat  over  the  interior  4  kHz  of  its 
passband.  A  couple  of  computer  routines  that  meet  these  requirements  admirably  are 
MATLAB’s  firls  and  firpm  routines,  where  the  Is  stands  for  least  squares  and  the  pm  stands  for 
Parks-McClellan.  A  4-kHz  passband  with  1-kHz  transition  bands  at  each  edge  is  easily  achieved 
by  either  of  these  routines. 

For  a  linear  filter  with  causal  impulse  response  duration  L  seconds,  the  transient  response  is 
of  length  L  seconds;  this  startup  transient  at  the  filter  output  must  be  discarded,  as  it  cannot  be 
fitted  properly.  In  addition,  the  time  delay  of  the  filtered  output  z,{nA)  is  LI 2  seconds,  not  L 
seconds;  accordingly,  excitation  x(n A)  must  be  delayed  by  this  same  LI 2  amount  before 
attempting  any  fits. 
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4.  SECOND-ORDER  CONSIDERATIONS 


4.1  UNIQUENESS  REQUIREMENT 


The  second-order  Volterra  component  waveform  was  given  in  equation  (1 )  as 


A  I  A'  I 


V2 (")  =  X  2 (*„ *2 )  x( n  -  *■ )  x{n -k2). 


(24) 


*1=0  *,=0 


It  must  be  noted  immediately  that  there  can  be  no  unique  solution  for  the  second-order  time- 
domain  kernel  h2(kl,k2).  Only  the  sum 


h2  (kt,k2)  +  h2(k2,kt) 


(25) 


of  symmetrically  located  kernel  values  (about  the  45°  line  in  A,,  A,  space)  affects  output  y2(n). 
This  is  due  to  the  symmetry  of  the  product  of  the  two  .v  values  in  variables  A', ,  As  in  equation 


(24).  To  enable  a  unique  solution  for  the  second-order  time-domain  kernel,  it  is  necessary  to 
impose  some  condition.  The  condition  adopted  here  is  to  take  the  second-order  kernel  to  be  real 
and  symmetric: 


h2(k2,kt)  =  h2  (A,,  A,). 


(26) 


There  is  no  loss  of  generality  in  imposing  this  symmetry  relation  on  the  second-order  time- 
domain  kernel  A, (A,,A2),  because  the  output  y2(n)  of  equation  (24)  depends  only  on  the  sum  in 
equation  (25),  no  matter  what  the  excitation  .v(/?)  is.  Also,  significant  advantage  will  be  taken  of 
this  symmetry  property,  including  a  reduction  in  the  number  of  unknown  coefficients  that  must 
be  determined,  as  well  as  avoidance  of  singular  matrices. 

4.2  PROPERTIES  OF  THE  SECOND-ORDER  FREQUENCY-DOMAIN  KERNEL 

The  second-order  frequency-domain  (complex)  kernel  corresponding  to  h2(klA,k2A)  is 


A'  I  A’  1 


(27) 


A.  -0  A--0 


The  inverse  relation  is 


(28) 
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Let  the  time  duration  of  h2(k{A,k2A)  in  each  time  dimension  be  T seconds,  and  let  the  frequency 
extent  of  //2(/,,/2)  in  each  frequency  dimension  be  (~F,F)  Hz.  These  limits  are  identical  to 
those  used  for  the  first-order  situation.  Then,  using  symmetry,  approximately  (2TF)2 / 2  samples 
are  required  to  completely  characterize  the  real  second-order  time-domain  kernel  h2(k]A,k2A). 
This  rapid  increase  in  the  number  of  required  coefficients  for  characterization,  in  progressing 
from  first  order  to  second  order,  is  the  COD. 

Since  h2(k]A,k2A)  is  forced  to  be  real,  it  follows  immediately  from  equation  (27)  that 


(29) 


And,  from  equations  (26)  and  (27),  there  follows 


h1U'2J\)  =  h2U\,J2). 


(30) 


That  is,  the  complex  second-order  frequency-domain  kernel  //2(/,,/2)  satisfies  both  a 
symmetry  relation  about  the  +45°  line  in  space,  and  a  conjugate  symmetry  relation 

through  the  origin  of  f^f2  space. 

These  symmetry  properties  mean  that  complex  frequency-domain  kernel  //,(  /,,/, )  needs 
to  be  specified  only  in  a  90°  sector  of  fvf2  space;  it  is  then  automatically  fixed  in  the  rest  of  the 
two-dimensional  frequency  plane.  In  particular,  the  “east”  sector,  composed  of  the  90°  sector 
centered  on  the  positive  f{  axis,  is  adopted  as  the  “fundamental  region”  in  ft,f2  space. 
Mathematically,  that  is  the  region 

f\  —  0,  -/,</,</,.  (31) 

Advantage  must  always  be  taken  of  these  symmetry  properties  because  they  afford  a  significant 
reduction  in  the  number  of  unknowns  that  must  be  determined  (estimated).  Only  the  values  of 
H2(f\,f2)  in  the  fundamental  region  (31)  need  to  be  specified. 

The  area  of  the  complete  frequency  space  /,,/2,  of  extent  (~F,F)  in  each  dimension,  is 
(2  F)2 .  On  the  other  hand,  the  area  covered  by  equation  (3 1 )  is  only 

\df]  )  df2=F2=~(2F)2.  (32) 

0  -/,  1  2- 

The  first  factor  of  1/2  is  due  to  use  of  the  conjugate  symmetry  relation  (29),  while  the  second 
factor  of  1/2  is  due  to  the  symmetry  relation  of  equation  (30).  These  factors  constitute  a 
significant  reduction  in  the  number  of  complex  coefficients  that  need  to  be  determined  at  second 
order. 
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A  seemingly  “obvious”  partitioning  of  /,  ,/2  spaee  is  to  take  squares  of  size  W  x  W.  If 
allowed,  this  would  reduce  the  COD  by  the  tremendous  factor  of  (F /W)1 .  However,  it  is  found 
that  the  outputs  of  these  individual  frequency  squares  interact  and  contaminate  each  other.  In 
addition,  it  is  not  obvious  how  the  off-diagonal  squares  in  which  cover  different 

frequency  ranges  in  /,  and  f2 ,  should  be  fit  to  a  single  band  of  frequencies  of  filtered  output 
z(nA). 


4.3  ALTERNATIVE  FORM  FOR  SECOND-ORDER  OUTPUT  y2(nA) 

Substitution  of  equation  (28)  into  equation  (24),  and  interchange  of  double  summation  and 
integration,  leads  to  the  relation 

v2(«A)  =  A2  [f  dfx  dj\  exp[/2/r(/,  +./>A]  //2(/,,/2)  X(f,)  X(f2).  (33) 

This  is  not  a  double  Fourier  transform;  there  is  only  one  time  variable  on  the  right-hand  side, 
namely,  «A. 

The  key  observ  ation  to  make  at  this  juncture  is  that  the  only  place  that  time  variable  ;;A 
appears  on  the  right-hand  side  of  equation  (33)  is  with  the  frequency  combination  J\  +  /, .  If 

second-order  Volterra  output  V2(nA)  is  to  have  frequency  eontent  only  in  the  band  for 

purposes  of  fitting  to  a  corresponding  filtered  version  of  z(nA),  and  if  X{  f)  is  broadband,  then 
second-order  frequency-domain  kernel  //2 (/,,/,)  must  be  restricted  to  be  nonzero  only  for 


f„  </,  +/:</,  (34) 

(and  the  corresponding  negative  frequencies).  This  condition  allows  the  complex  exponential  in 
equation  (33)  to  take  on  frequency  variation  only  in  the  band  (fa,fh).  The  region  in  equation 

(34)  is  definitely  not  square  in  space.  Rather,  see  the  blue  and  green  regions  in  figure  2. 

Equation  (34)  describes  an  infinite  strip  at  angle  -45°  in  the  plane,  with  perpendicular 
width  (fh  -  /u)/V 2  =  Wj  yfl .  However,  the  fundamental  region  is  limited  to  be  below  the  +45° 
line  in  the  /, ,/,  plane;  see  the  red-bordered  90°  east  region  in  figure  2.  In  addition,  frequency 
fx  cannot  exceed  the  limit  F.  The  shape  of  this  finite  confined  strip  in  the  /,,/,  plane  is  similar 
to  the  shape  of  the  state  of  Nevada.  This  is  the  restricted  region  of  /,  ,/2  spaee  in  which 
H2(f„f2)  is  allowed  to  be  nonzero  if  y2(nA)  in  equation  (33)  is  to  contain  frequency  content 
limited  to  the  frequency  range  (/„,//,). 
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Figure  2.  Diagonal  Strips  in  fij~2  Plane 


Frequency  fa  can  start  at  zero.  Also,  it  is  not  necessary  to  let  frequency  fb  exceed  F.  It  is 
assumed  that  there  is  no  frequency  content  in  nonlinearity  output  z(nA)  beyond  F.  Thus,  the 
fundamental  (blue)  strip  can  slide  anywhere  between  the  dashed  red  and  black  lines.  The  four 
yellow  Xs  represent  the  locations  of  the  symmetry  points  of  //,(/ , ,/,)  dictated  by  equations 
(29)  and  (30). 

The  length  of  the  frequency  strip  (34)  along  the  -45°  line  starts  at  value  ~Jl  F  for  fa  =  0. 

Therefore,  the  area  covered  by  each  strip  is  approximately  -Jl  F  x  w/Jl  «  F  W,  not  F2  as  it 
would  be  if  frequency  partitioning  were  not  employed  at  all.  Thus,  the  COD  has  been  alleviated 
(not  eliminated)  by  a  factor  of  W/F  through  frequency  partitioning,  a  very  worthwhile  and 
welcome  reduction  in  many  cases.  Since  the  frequency  “area”  covered  by  each  elemental 
complex  exponential  in  equation  (17)  is  (1  jT)~ ,  the  number  of  unknown  coefficients  is  about 
FW j (\/T)2  =  TF  TW,  not  {TF)1 . 
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4.4  ALTERNATIVE  DERIVATION  OF  RESULT  (34) 

A  Fourier  transformation  of  equation  (33)  results  in  spectrum 

n(/)=  J#i  tf2(/r/-/>) *(/,)*(/-/,)•  (35) 

In  order  for  Y2(f)  to  he  nonzero  only  in  frequency  band  (/„,/* ),  function  H2(J\,f-f\)  must 
be  nonzero  only  for  fa  <  f  <  fh.  Let  f2  -  f  -  fv  Then,  /  =  J\  +  /,,  from  which  it  then 
follows  that  condition  (34)  must  be  satisfied. 


4.5  BASIS  FUNCTIONS  IN  TWO  DIMENSIONS 


The  basis  functions  in  one  dimension  (first  order)  were  derived  in  equations  (17)  and  (18), 
starting  from  an  elemental  complex  exponential  in  the  fundamental  range.  Here,  in  two 
dimensions  (second  order),  the  corresponding  complex  elemental  starting  function  is 


[a(//J, ,  m2 )  -  i  b(iiit ,  m2 )]  exp  ;  2  n  —  A',  A  + i  2  n  ■— -  A,  A 


(36) 


where  coefficients  {a(mt ,/«,)}  and  {b(mt,m2)\  are  real,  and  frequencies  fx  =  mjT  and 
f2  =  m2/T  are  confined  to  the  fundamental  (blue)  region  in  the  /, ,f2  plane,  as  described 
earlier.  The  total  real  function  for  expanding  second-order  time-domain  kernel  /?,  (A,  A,  A,  A) , 
using  the  conjugate  and  reflective  symmetries  in  equations  (29)  and  (30),  respectively,  is  then 
(with  scaling  1/4  =  1/2  *  1/2!) 


—[a(m^m2)-i  )]  exp(i  A)  +— [a(mt,w2 )  —  i b(m^, ///,)]  exp(/  B ) 

+  complex  conjugate  terms 

(37) 


=  «(«.,/»,) 


cos(A)  +  cos(B) 


+  b(m.,m2) 


sin(/l)  +  sin(£) 


where 


A  = 


In  Is. 
T 


(//»,  A,  +  m2  A; ), 


B  = 


2n  A 
T 


( 117 2  A |  +  /77|  A2  ) . 


(38) 


The  two  basis  (bracketed)  functions  in  the  second  line  of  equation  (37)  are  obviously  real  and 
symmetric  in  A,,A2 ,  upon  use  of  equation  (38),  as  expected  and  required.  However,  their  forms 
are  somewhat  surprising  and  not  easy  to  anticipate,  without  the  construction  utilized  in  equation 
(36)  and  the  first  two  lines  of  equation  (37). 
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It  was  pointed  out  in  the  sequel  to  equation  (15)  that  the  coefficient  b( 0)  in  the  first-order 
case  was  absent  because  the  corresponding  sine  term  in  the  expansion  was  zero  for  m  =  0.  A 
more  subtle  but  related  property  holds  at  second  order.  From  equations  (37)  and  (38),  a  typical 
term  in  the  sine  basis  set  is 


—  sin 
2 


(2nA,  /  /  ^ 

- (m,  k.  +  m1  k-, ) 

T 


+  — sin 
2 


2n  A  N 

— — (w2  A:,  +m]k2) 


V 


y 


(39) 


This  term  is  obviously  zero  for  the  point  aw,  =  0,  m2  =  0  in  fx,f2  space.  But  equation  (39)  is 
also  zero  for  /»,  +  m2  =0,  even  though  neither  term  is  zero.  This  latter  condition  yields  a  line  in 
the  two-dimensional  frequency  plane  All  of  the  terms  corresponding  to  this  line, 

aw,  +  m2  =  0,  must  be  eliminated  from  the  sine  basis  set  in  order  to  avoid  a  singular  matrix  in  a 
second-order  fitting  process.  None  of  the  [cos(/l)  +  cos(£)]/2  terms  in  equation  (37)  is  zero. 


4.6  SECOND-ORDER  KERNEL  EXPANSION 

Combine  equations  (37)  and  (38)  to  form  the  fundamental  second-order  basis  functions: 


1 

—  cos 

ln  i  i- 
- (aw,  k. 

+  tn  j  ) 

1 

+  —  cos 

- (//?.  k-,  +  k, ) 

2 

\kk  1  1 

2 

K  1  * 

(40) 


s2(rnl,m2;ki,k2)  =  —  sin 


— —  (///,  k^  +  m2  k2) 
K 


1  . 

+  —  sin 
2 


k:  +  m2  k\ ) 


using  T  =  KA.  Both  of  these  basis  functions  are  symmetric  in  kx,k2,  and  are  symmetric  in 
mx ,  m2 .  Also, 


c2(0,0;A:i  ,k2)  *  0;  s2(mx,-m[;kl,k2)  =  0  for  all  kx,k2. 


(41) 


bor  notational  convenience,  define 


2  n 

2n 

COS 

- ink 

,  s(n/,k)  =  sin 

- ink 

K  J 

U  J 

(42) 
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Then,  equation  (40)  can  be  expanded  as 


c2(mx,m2\kx,k2)  =  —  [c(mx,kx)  c(m2,k2)-s(mx,kx)  s(m2,k2) 
+  c(aaai,A2)c(aw,,A1)-s(aaa,,A2).s(aa*2,Ai)], 

s2(mx,m2‘,kx,k2)  =  [s(aaa, , kx )  c(m2,k2)  +  c(mx,kx )  s{m2,k2 ) 
+ ,?(/«, , k2 )  c(m2 , kx )  +  c(/w, ,k2)  s(m2 , A, )]. 


(43) 


Now,  referring  back  to  equation  (37),  the  expansion  for  the  second-order  kernel  is  given  by 


h2(kx\k2A)=-YJ  a(mxtm2)  [c(mx,kx)  c(m2,k2)-s(mx,kx)  s(m2,k2) 

2  Mc 

+  c(rnx , k2 )  c(m2 ,  A, )  -  s(mx ,  A, )  s(m2 ,  A, )] 

+  z:  Yj  ^('"i  'mi )  [■'(/«, , A, )  c(m2 , A, )  +  c(aaa,  , A, )  s(m2,k2) 

+  s(mx  ,A2 )  c(a?/2  ,A, )  +  c(a//,  ,  A, )  ,v(a7/2  ,A, )], 

where  A/(  and  Ms  correspond  to  the  particular  two-dimensional  region  of  frequency  space 
fx,f2  in  figure  2  that  is  under  investigation.  Region  Ms  is  smaller  than  Mc  by  virtue  of  the 
comments  following  equation  (39).  Second-order  kernel  /a2(A,A,A2A)  in  equation  (44)  is 
symmetric  in  A,, A,. 


4.7  SECOND-ORDER  WAVEFORM  EXPANSION 

As  in  equations  (22)  and  (42),  define  sequences 

K  I 

x,  (n,m)  =  Y  x(n-k)  c(m,k), 

(45) 

K  1 

xs  (a/.aaa)  =  Y  x(n  -  k )  s(m,k). 

*=0 

Sequence  {^(aj.aw)}  is  even  in  aw,  while  {x  (//,//?)}  is  odd  in /?/.  Also,  .vs(//,0)  =  0  for  all  n. 

These  convolution  sequences  can  be  calculated  just  once  and  stored  in  two  N  X  M  arrays  for 
further  processing  as  needed.  Substitution  of  equation  (44)  into  equation  (24),  interchange  of  the 
double  summations,  and  the  use  of  equation  (45)  yields  the  second-order  Volterra  component  in 
the  form 
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(46) 


y  2  (  n)  =  Z  ’  n,2  )  t*r  '"l  )*<-(«’  ,,J2  )~Xs(n’  m\  )*,(">  »>2  )1 

Mc 

+  X  ^(/w,  i  ffli )  K  («,"*,)  (/?,  /«, )  +  xc  (n, mx )  xs  (//.  m2 )] , 

ws 

where  duplicate  terms  have  been  added  together.  Both  bracketed  terms  are  symmetric  in  /?/,  ,///2. 
This  last  equation  for  y2(n)  is  now  in  a  form  ready  for  use  in  a  least  squares  fit  to  a  filtered 
version  of  nonlinearity  output  z(n). 


4.8  DISCARDING  OF  EDGE  STRIPS  AT  SECOND  ORDER 

In  an  earlier  first-order  subsection,  a  method  for  the  discarding  of  edge  coefficients  for  the 
first-order  fitting  procedure  was  presented.  That  method  can  be  generalized  to  second  order  as 
follows.  The  1  -kHz  diagonal  strip  next  to  the  bottom  left  edge  of  each  6-kHz  diagonal  strip  in 
the/1,/2  plane  must  be  discarded.  That  particular  1-kHz  diagonal  strip  in/,/  leads  to  the 
frequency  components  of  y2{nA)  at  the  lower  edge  of  the  (/,/>)  band,  namely,  to  components  in 
the  band  (fa,fa  +  1  kHz). 

Similarly,  the  1-kHz  diagonal  strip  next  to  the  upper  right  edge  of  each  6-kHz  diagonal  strip 
must  be  discarded.  Only  the  intenor  4-kHz  diagonal  strip  of  each  6-kHz  diagonal  strip  is 
retained.  When  the  adjacent  frequency  band  comes  under  investigation  for  fitting  purposes, 
overlap  will  again  be  required.  But  the  next  analysis  strip  will  be  moved  over  by  4  kHz,  not  6 
kHz.  Thus,  two  edge  bands  are  discarded  at  first  order,  but  two  edge  strips  are  discarded  at 
second  order.  Third  order  will  likely  require  discarding  of  edge  volumes  in  the  three- 
dimensional  frequency  space/,/,/. 


4.9  TEST  PROCEDURE  FOR  A  FIRST-ORDER  AND  SECOND-ORDER  EXAMPLE 

A  nonlinear  system  with  known  first-order  kernel  hXe(k\  A)  and  known  second-order  kernel 
lt2e(klA,k2A)  was  excited  with  stationary  random  process  x(nA).  The  output  z,  (//A)  of  the 
first-order  kernel  hu  was  computed  and  stored.  Also,  the  output  z,(//A)  of  the  second-order 
kernel  h2e  was  computed  and  stored.  Finally,  the  total  nonlinear  system  output 

z(»A)  =  z0 +z,(//A)  +  z2(»A)  (47) 

was  computed  and  stored,  as  was  colored  excitation  x(nA).  See  figure  3. 
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NONLINEAR  SYSTEM 


Figure  3.  Test  Procedure  for  Second  Order 


Working  only  with  the  recordings  of  excitation  x(nA)  and  total  nonlinear  response  z(/?A), 
the  object  was  to  choose  real  coefficients  a(m)  and  b{m)  for  first-order  Volterra  kernel  /?,  (A:,  A), 
and  to  simultaneously  choose  real  coefficients  a(mx,m2)  and  b(m^,m2)  for  second-order 
Volterra  kernel  //,(£,  A,  A: A),  such  that  the  total  Volterra  model  output 

v(  nA)  =  vn  +  v’|  (nA)  +  y2  ( nA)  (48) 

fitted  the  jittered  waveform  zab(nA)  as  well  as  possible  in  a  least  squares  sense. 

Five  checks  can  be  conducted  on  the  estimates  obtained  from  this  control  example.  Namely, 
one  can  compare:  estimate  /;,  with  the  exact  hie;  estimate  h2  with  the  exact  h2e\  estimate  y\ 

with  the  exact  z,;  estimate  v2  with  the  exact  z2;  and  estimate y  with  the  exact  z.  In  an  actual 
practical  situation  of  an  unknown  nonlinear  system  under  investigation,  the  only  available 
comparison  will  be  least  squares  model  response y  with  nonlinearity  response  z.  This 
comparison  can  also  be  conducted  band-by-band  in  frequency. 
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4.10  INTERPRETATION  AND  USEFULNESS 


The  power  in  the  total  model  output  j(«A)  relative  to  the  power  in  total  nonlinear  output 
z(nA)  measures  the  goodness  of  the  fit.  If  only  a  small  fraction  of  the  z(/7A)  power  is  accounted 
for,  there  may  be  independent  noise  generated  inside  the  nonlinear  system  that  can’t  possibly  be 
fit;  the  level  of  this  internal  noise  can  be  determined  by  turning  off  input  excitation  x(«A)  and 
observing  the  inherent  power  in  output  z(nA).  Or  the  order  of  the  Volterra  (or  Wiener)  model 
may  simply  be  inadequate  and  must  be  increased  for  a  better  fit.  This  would  indicate  that  higher 
order  nonlinearities  are  present  in  the  system  under  investigation. 

The  relative  levels  of  the  powers  in  the  total  final  fits  (nA)  and  y2(nA)  can  be  used  as 
indicators  of  the  relative  amount  of  second-order  nonlinearity  in  the  system  of  interest.  If  the 
nonlinear  power  is  significant  in  comparison  with  the  linear  power,  it  will  behoove  one  to 
account  for  this  nonlinearity  in  the  system  in  the  future. 

Volterra  kernel  estimates  /7,(A,A)  and  /72(A,A,A2A)  can  now  be  used  to  predict  the 
performance  of  the  investigated  nonlinear  system  to  different  future  inputs,  without  actually 
having  to  perform  those  experiments. 
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5.  SECOND-ORDER  NUMERICAL  RESULTS 
FOR  THE  CONTROL  EXAMPLE 


The  following  results  were  accomplished  by  frequency  partitioning.  Figure  4  displays  the 
results  for  the  first-order  kernel.  The  exact  kernel  hu,{kx )  =  A,,  (A, )  is  drawn  in  blue  (b),  and  the 

estimated  model  kernel  A,  (A,)  is  drawn  in  black  (k).  The  error  eAI(A,)  between  them  is  drawn 

in  red  (r).  The  estimate  is  a  virtual  ov  erlay  of  the  exact  kernel;  there  is  a  small  region  at  the  far 
right  w  here  the  red  error  curve  deviates  slightly  from  zero.  The  ratio  of  the  energy  in  the  error 
sequence  to  the  energy  in  the  exact  kernel  sequence  is  0.000042. 

The  corresponding  results  for  the  second-order  kernel  h2i, ( A, , A2 )  =  h2: (A,  ,k2)  are  given  in 
figures  5  and  6.  The  estimated  model  kernel  in  figure  6  is  visually  indistinguishable  from  the 
exact  kernel  in  figure  5.  The  ratio  of  the  energy  in  the  error  between  the  two  functions  to  the 
energy  in  the  exact  kernel  is  0.00023. 

Figure  7  gives  the  results  for  a  section  of  the  first-order  waveforms.  The  error  w  aveform 
ex(n)  between  estimate  yt (n)  and  exact  waveform  z,(/?)  is  extremely  small  and  has  a  relative 
error  energy  of  0.000043,  the  same  as  for  the  first-order  kernels  in  figure  4. 

The  corresponding  comparison  for  the  second-order  waveforms  is  presented  in  figure  8. 

The  variation  with  time  of  these  second-order  waveforms  is  much  more  erratic  and  quicker  than 
for  the  first-order  results  in  figure  7.  This  is  not  unexpected  for  a  second-order  nonlinearity.  The 
ratio  of  error  energy  to  exact  energy  is  0.000 1 4. 

Finally,  the  total  waveform  results  for  estimate  y(n)  and  nonlinearity  response  waveform 
z{n)  arc  given  in  figure  9.  The  error  energy  ratio  is  0.000042  for  this  case.  This  last  plot  is  the 
only  comparison  that  can  be  constructed  for  the  practical  examples  where  only  input  excitation 
x(n)  and  nonlinear  response  z(n )  are  available. 

Two  new  concepts  have  been  introduced  at  second  order  in  this  investigation.  They  arc  the 
fundamental  region  (31)  and  the  special  basis  functions  (40)  in  two-dimensional  time  space 
A|,  A'2-  These  items  are  crucial  in  achieving  a  satisfactory  second-order  fit  by  means  of  least 
squares.  Also,  it  was  then  necessary  to  eliminate  any  zero  basis  functions,  and  to  use  overlap  and 
to  discard  edge  estimates  that  were  untrustworthy.  The  combination  of  all  these  items  then  led  to 
a  viable  and  accurate  method  of  estimating  the  internal  behavior  of  a  second-order  nonlinear 
system  with  memory. 
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amplitude 


Figure  4.  First-Order  Kernels  hji(b)  hi(k)  e/,i(r) 


400  400 


*1 

Figure  5.  Exact  Second-Order  Kernel  li22(k],kx ) 
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Figure  6.  Model  Second-Order  Kernel  /i2(k/,kz) 


Figure  7.  First-Order  Waveforms  zi(b)  }'i(k)  efr) 
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amplitude  amplitude 


Figure  8.  Second-Order  Waveforms  12(b)  }'i(k)  ei(r) 


Figure  9.  Waveforms  z(b)  y(k)  e(r) 
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6.  THIRD-ORDER  CONSIDERATIONS 


6.1  UNIQUENESS  REQUIREMENT 

The  third-order  Volterra  component  waveform  was  given  in  equation  (1 )  as 

=  2  X  Z  fh(knk2’k3)  x{n-kx)  x(n-k2)  x(n-k2).  (49) 

*,=0  *3=0  *,=0 


It  must  be  noted  immediately  that  there  can  be  no  unique  solution  for  the  third-order  time- 
domain  kernel  /j3(A,,A2,A3).  Only  the 


/2j(A|  ,  A'2 , ^3 )  +  /t3 (Ar3 ,  A | , A2 )  +  hy (A2 , A3 , A( )  +  h2 (ki , A3 , A2 )  +  h2 (k2 , Aj ,  A3 )  +  /j3 ( A, , A2 , k\ )  (50) 


of  six  symmetrically  located  kernel  values  in  A,  ,A2,A3  space  affects  output  v3(/j).  This  is  due  to 
the  symmetry  of  the  product  of  the  three  x  values  in  variables  A,,A2,A3  in  equation  (49).  To 

enable  a  unique  solution  for  the  third-order  time-domain  kernel,  it  is  necessary  to  impose  some 
condition.  The  condition  adopted  here  is  to  take  the  third-order  kernel  to  be  real  and  symmetric: 


h} (A| , A2 , A3 )  —  /?3(A3,A|,A2)  —  /j3 ( A2 , A3 , A| )  —  /t3(A,,A3,A2)  —  /?3 {k2 , A| , A3 )  —  h2 (A3, A2 , Aj ). (5 1 ) 


There  is  no  loss  of  generality  in  imposing  this  symmetry  relation  on  the  third-order  time- 
domain  kernel  hx(k{  ,A2,A3),  because  the  output  y}(n)  of  equation  (49)  depends  only  on  the  sum 
in  equation  (50),  no  matter  what  the  excitation  .y(»)  is.  Also,  significant  advantage  will  be  taken 
of  this  symmetry  property,  including  a  reduction  in  the  number  of  unknown  coefficients  that 
must  be  determined,  as  well  as  avoidance  of  singular  matrices. 


6.2  PROPERTIES  OF  THE  THIRD-ORDER  FREQUENCY-DOMAIN  KERNEL 

The  third-order  frequency-domain  (complex)  kernel  corresponding  to  /j3(A|A,A2A,A3A)  is 


x  i  x  i  a:  i 


//,(/,,/,,/,)=  X  Z  Z  ^i(^iA,A,A,A3A)  exp(-/2;r  /,  A, A  -  i2n  f2  A2A  -  iln  f2  A3A).  (52) 

*  =0  *1=0  *,=() 

The  inverse  relation  is 

/i3(A,A,A2A,A3A)  =  A3  JJJV/’i  df2  df2  exp(/2^/,  A,A  +  /2;r  f2  k2A  +  i2n  J\  A3A).(53) 
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Let  the  time  duration  of  h3(ktA,k2A,k3A)  in  each  time  dimension  be  T seconds.  And  let  the 
frequency  extent  of  H3{ /,,/2,/3)  in  each  frequency  dimension  be  (—F,F)  Hz.  These  limits 
are  identical  to  those  used  for  the  first-  and  second-order  situations.  Then,  using  symmetry, 
approximately  (2TF)X  /b  samples  are  required  to  completely  characterize  the  third-order  time- 
domain  kernel  h-li(k]A,k2A,k2iA).  This  rapid  increase  in  the  number  of  required  coefficients  for 
characterization,  in  progressing  from  first  order  to  second  order  to  third  order,  is  the  COD. 

Since  hi(kiA,k2A,k3A)  is  forced  to  be  real,  it  follows  immediately  from  equation  (52)  that 

~A~A)  -  «,(/„/,./,)*•  (54) 


And,  from  equations  (51)  and  (52),  there  follows 


=  #  3  (/, .  /, ,  f2 ) =  H3  (f2  ,/„/,)  =  //,  (/3 ,  f2 ,  fx ). 


(55) 


Thus,  the  complex  third-order  frequency-domain  kernel  //,(/,,  /2,  /3)  satisfies  a  conjugate 
symmetry'  relation  through  the  origin  of  three-dimensional  /j,/2,/3  space,  as  well  as  numerous 
symmetry  relations. 


6.3  ALTERNATIVE  FORM  FOR  THIRD-ORDER  OUTPUT  y3(nA) 

Substitution  of  equation  (53)  into  equation  (49),  and  interchange  of  triple  summation  and 
integration,  leads  to  the  relation 

y3(n A)  =  A1  jj{#,  df2  JJ\  exp[i2x(f,  +/2  +/,)i»A] H3(ft,f2,f3)  X{j\)  X(f2)  X{f3).{ 56) 

This  is  not  a  triple  Fourier  transfonu;  there  is  only  one  time  variable  on  the  right-hand  side, 
namely,  nA. 

The  key  observ  ation  to  make  at  this  juncture  is  that  the  only  place  that  time  variable  nA 
appears  on  the  right-hand  side  of  equation  (56)  is  with  the  frequency  combination  /,  +/,  +/3. 

If  third-order  Volterra  output  y3(nA)  is  to  have  frequency  content  only  in  the  band  (/„,/<,),  for 
purposes  of  fitting  to  a  corresponding  filtered  version  of  z{nA),  and  if  X(f)  is  broadband,  then 
frequency-domain  kernel  //3(/,,/2,/3)  must  be  restricted  to  be  nonzero  only  for 


.fa  </l  +fl  +f3  <fh 


(57) 
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(and  the  corresponding  negative  frequencies).  This  condition  allows  the  complex  exponential  in 
equation  (56)  to  take  on  frequency  variations  only  in  the  band  ).  The  lower  bound  in 
equation  (57)  describes  a  plane  in  three  dimensions,  while  the  upper  bound  describes  a  parallel 
plane.  The  in-between  region  in  equation  (57)  is  an  infinite  tilted  slab  in  three  dimensions. 

The  restricted  frequency  region  at  second  order  was  given  in  equation  (3 1 ).  That  result  can 
be  extended  to  third  order  by  observing  that  only  one  of  the  six  equal  third-order  frequency- 
domain  values  in  equation  (55)  needs  to  be  specified.  In  particular,  one  selects  as  the  primary 
ease  that  which  satisfies  the  rule 

|/>N/;M/,|-  (58) 


This  limited  three-dimensional  frequency  region  is  further  reduced  by  application  of  equation 
(54),  leading  to  the  final  restricted  region 

0  <f<F.  (59) 

If  frequency  f\  is  held  fixed  at  a  positive  value,  this  allowed  region  in  the./:,/?  plane  consists  of 
two  90°  sectors,  namely,  the  east  and  west  sectors.  Thus,  the  shape  of  the  restricted  volume  in 
fufiifc  space  is  a  pair  of  triangular  pyramids  with  one  common  edge. 

Advantage  must  always  be  taken  of  these  symmetry  properties  because  they  afford  a 
significant  reduction  in  the  number  of  unknowns  that  must  be  determined  (estimated).  Only  the 
values  of  in  the  limited  region  (59)  need  to  be  specified. 

The  volume  of  the  complete  frequency  spaee/1,/2,/3,  of  extent  (~F,F)  in  each  dimension, 
is  (2 F)\  On  the  other  hand,  the  volume  covered  by  equation  (59)  is  only 

)  df2  }  df2  2|/2|  =  J J/i  If;  =  =  i  i (2 F)\  (60) 

The  factor  of  1/2  is  due  to  use  of  the  conjugate  symmetry  relation  (54),  while  the  factor  of  1/6  is 
due  to  the  symmetry  relations  of  equation  (55).  These  factors  constitute  a  significant  reduction 
in  the  number  of  complex  coefficients  that  need  to  be  determined  at  third  order. 


6.4  BASIS  FUNCTIONS  IN  THREE  DIMENSIONS 


Corresponding  to  the  two-dimensional  development  in  equation  (36),  the  complex  elemental 
starting  function  in  three  dimensions  (third  order)  is 


[a(m,,rn2,m} )  - 1  />(//;, ,///,,/»,)]  expj  i27T^-k,A  +  i 2n^j^k2A  +  i2n'-^rkyA 


(61) 
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where  coefficients  {a(mi, m2,mi)}  and  {b(m[,m2,mi)}  are  real,  and  frequencies 


fi  =m\ /T*  fi  =mi/T’  fi  =mi/T 


(62) 


are  confined  to  the  region  in  the  /,,/2,/3  plane  described  in  equation  (59).  Let 


a  =  2x A/T  =  2x/K,  c  =  a(ml,m2,mi)  —  /  b(mx,m2,m2). 


(63) 


Then,  using  equation  (61 ),  the  total  elemental  real  function  for  expanding  third-order  time- 
domain  kernel  /7,(A,A,A2A,  A3A),  using  the  conjugate  and  reflective  symmetries  in  equations  (54) 

and  (55),  respectively,  is,  after  expansion,  simplification,  and  scaling  by  1/12  =  1/2  *  1/3!, 

-|V  { c  exp[/a  (/n,Ar,  +m2k2  +  miki)]  +  c  exp[/6r  (/n3A,  +  m[k2  +m2k^)] 

+  c  exp[/  a  (m2k{  +  mik2  +  mlki )]  +  c  exp[/  a  ( mlk]  +  myk2  +  tn2ki )] 

+  c  exp[/«  (ni-,kl  +mxk-,  +  miki)]  +  c  exp[/«  (mikl  +m-,k 2  +  7W|A3)]} 

-  -  ...  (64) 

+  complex  conjugate  terms 


=  a(m{ ,  m2 ,  my )  va  +  b(mx ,  m2 ,  m3 )  vh , 
where  the  real  basis  functions  are 

va  =j  {cos[«  (mlkl  +m2k2  +??23A-3)]  +  cos[a(m3/:l  +  mxk2  +  ;m,A3)] 

+  cos[a  (///,£,  +  myk2  +  m,/:3)]-i-cos[a(ml/:l  +tnik2  +  ?w2A3)] 

+  cos[a(???;A,  +  mxk2  +?773  A3)]  +  cos[a(/?2,A'l  +  m2k2  +  ;?;,A3)]}, 

vb  =  •£  {sin[a(A7/,Al  +m2k2  +/w3A3)]  +  sin[«(/H,Al  +mxk2  +  /w2A3)] 

+  sin[«(??22Al  +m3k2  +  ?wlA3)]  +  sin[Q'(??71A,  +m2k2  +  ??/,A3)] 

+  sin[«(?w,Al  +  mlk2  +  /w3A3)]  +  sin[#(/n3Al  +m2k2  +/m,A3)]}. 

To  minimize  computer  storage  requirements,  both  functions  in  equation  (65)  must  be 
expanded  and  expressed  in  terms  of  the  quantities  c(m,k)  and  s(m,k)  defined  in  equation  (42). 
The  end  results  are 

v„  =i[c(ml,Ar,)/l  -s(ml,kl)t2  +  c(mi,k2)  /3  -s(m{  ,A2)  t4  +  c(mx,k})t5  -.v(;n,,A3) /J, 


(66) 
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where 


/,  =  c(/7?, ,  k2 )  c{my ,  A3 )  -  s(m2 ,  k2 )  5(77?, ,  A, )  +  c(m2 ,  A3 )  c(my ,  A, )  -  s(m2 ,  A, )  s(m3 ,  A, ) , 
t2  =  s{m2 ,  A, )  c(m3 , k3 )  +  c(w2 ,  k2 )  s(m3 , k3 )  +  s{  //?, , k3 )  c( m3 ,k2)  +  c(m2 , k3 )  s( m3 ,A,), 
t3  =  c(w2 , k3 )  c{m3 ,  k] )  - s( m2 ,  k3 )  5(777, , kx )  +  c(m2 ,  A, )  e(/w3 ,  A, )  -  5(77?, ,  A, )  s(m3 ,k3), 
tA  =  s(m2 ,  k3)  c(m3 ,  A, )  +  c(m2 ,  A, )  5(ot3,  A, )  +  5(777., ,  A, )  <:(//?, ,  A, ) + c(w, ,  A, )  a(»?3 ,  A, ) , 
t5  =  c(???2 ,  A, )  c(m, ,  A, )  -  5(7/?, ,  A, )  5  ( m3 ,  A, )  +  c(/77, ,  A; )  c(;w3 ,  A, )  -  .v (m2 ,  A, )  5(777, ,  A, ) , 
h  -  s( m2  >  *1 )  C(mi  7  Aj  )  +  c(w2 ,  A,  )  i'(/773 ,  A;  )  +  5(777,  ,  A,  )  c(/77, ,  A, )  +  C'(777,  ,  A,  )  5( 777, ,  A,  )  . 


6.5  THIRD-ORDER  KERNEL  EXPANSION 

Thus,  from  equation  (64),  the  expansion  for  the  third-order  kernel  is  given  by 
/73(A„A2,A3)  =  £  a(m,,m2,m3)  v,  +  ^  6(777,, 777,, 777,)  v6 ,  (68) 

.W,  M, 


where  vo  and  are  given  by  equations  (66)  and  (67).  Region  Mc  is  constituted  by  those 
//?, ,777,, 777,  values  in  equation  (62)  that  satisfy  the  conditions  in  equation  (59).  Region  A/(  is 
identical  to  Mc  except  that  the  points  /??,  =  0,  /??,  +  m3  =  0  and  m2  =  0, 777,  +  7 ??,  =  0  and 
777,  =  0, 777,  +  7/7,  =  0  must  be  discarded  because  these  particular  basis  functions  vh  in  equation 
(65)  are  7ero  for  all  A, ,  A, ,  A, . 


6.6  THIRD-ORDER  WAVEFORM  EXPANSION 


The  third-order  Volterra  component  waveform  is  now  given  by  substitution  of  equation  (68) 
into  equation  (49).  Af  ter  expansion,  simplification,  and  use  of  definitions  (42),  the  result  is 

V,  ( 77 )  =  ^  a ( ,  m2 , 777,  )  |>(  ( 77,  777,  )  C(77,  777,  ,  777,  )  -  X,  (77,  777,  )  5  (  77, 777,  ,  777, )] 

\f 

(69) 

+  ^  6(777, ,  7/7,  ,  777,  )  [*,  (77,777,  )  C(/7, 777,  ,777,  )  +  X( .  (/7,  777,  )  5(77,  777,  ,777, )] , 

M, 

where  sets  of  six  duplicate  terms  have  been  added  together,  and 


C(/7,//7,,//7,)  =  Xf  (77,  777,)  Xr  (77,  777,)  -  X,  (77,  777,)  X,  (77,  777,  ), 
5(77,  7/7,,  7/7,)  =  X,  (77,777,)  X,  (//,  7/7,  )  +  X,  ( 77 ,  777,  )  X,  (//,  777,  ). 


(70) 
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Both  bracketed  terms  in  equation  (69)  are  symmetric  in  and  involve  only  the  known 

stored  arrays  xc  and  xs .  Equation  (69)  is  now  in  a  form  ready  for  use  in  a  least  squares  fit  to  a 
filtered  version  of  nonlinearity  output  z(n). 
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7.  APPLICATION  TO  A  CICADA  SONG 


7.1  WAVEFORMS  AND  SPECTRA 

A  time  section  of  a  cicada  measurement  (number  SI  4)  is  displayed  in  figures  10  and  1 1 . 

The  data  record  is  about  4  seconds  long  and  the  sampling  frequency  is  fs  =  96  kHz.  The  plots  in 
figures  1 0  and  1 1  are  of  laser  measurements  L  \  and  Lz  of  the  tymbal  movements  on  both  sides  of 
the  cicada’s  body.  Figure  12  pertains  to  a  microphone  M  located  a  few  feet  from  the  insect.  The 
data  segment  to  be  analyzed  here  stretches  from  time  index  n  =  2.5e5  to  n  =  3.5e5,  which 
corresponds  to  approximately  1  second  of  data.  From  observ  ation  of  the  microphone  data  in 
figure  1 2,  it  is  obvious  that  this  is  a  noisy  measurement.  For  example,  although  the  two  lasers 
are  giving  zero  readings  for  n  =  0.4e5  to  n  =  1 ,3e5,  the  microphone  output  is  distinctly  nonzero 
and  noisy  during  this  time  interval.  Therefore,  it  is  expected  that  this  noise  continues  during  the 
1 -second  analysis  interval  indicated  above,  and  will  constitute  a  portion  of  output  data  z{n)  that 
cannot  be  fitted  by  a  Volterra  expansion.  The  laser  measurements  are  considered  as  two  inputs, 
X\(n)  and  xz(n),  to  a  nonlinear  system  with  memory  (the  cicada  body),  while  the  microphone 
measurement  is  considered  as  the  output  z(n)  of  that  same  nonlinear  system. 

A  blown-up  section  of  duration  0.05  second  of  laser  Lz  is  presented  in  figure  13.  Although 
very  noisy  and  jagged,  there  is  obviously  a  “periodic  behavior”  of  approximately  0.01  second.  A 
sample  time  segment  of  0.01  second  is  displayed  in  figure  14.  The  first-order  probability  density 
of  the  complete  laser  recording  is  approximately  exponential,  not  Gaussian.  Also,  some  very- 
high-frequency  components  are  present  in  the  input  data. 

The  spectra  of  lasers  L\  and  Lz  are  plotted  in  figures  15  and  16.  T  he  periodic  low-frequency 
component  shows  up  as  a  narrowband  component  centered  at  96  Hz,  as  expected.  Both  spectra 
have  a  great  deal  of  low-frcquency  content  and  are  fairly  flat  almost  all  the  way  out  to  the 
Nyquist  frequency  of  48  kHz.  It  appears  that  the  cicada  is  capable  of  generating  some  very-high- 
frequency  components,  all  the  way  up  to  50  kHz  or  so.  If  that  is  the  case,  the  96-kHz  sampling 
rate  may  not  have  been  large  enough  to  avoid  some  aliasing.  Recall  that  the  laser  time 
recordings  themselves  in  figures  10  and  1 1  are  very  clean;  that  is,  there  is  very  little  noise  in 
these  laser  recordings.  Thus,  one  has  a  non-Gaussian  and  non-white  process  to  use  for  excitation 
x(n).  This  situation  is  not  what  one  would  choose  for  an  excitation  if  the  input  were  under  one’s 
control. 

The  spectrum  for  the  microphone  recording  is  given  in  figure  1 7.  There  are  some 
pronounced  humps  at  6  kHz  and  8  kHz  that  bear  investigation.  There  is  a  strong  (unexplained) 
very  narrowband  component  centered  at  39.472  kHz.  And,  again,  the  spectrum  remains  flat 
almost  out  to  the  Nyquist  frequency  of  48  kHz. 
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Figure  10.  Waveform  Li 
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Figure  11.  Waveform  L2 
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Figure  12.  Waveform  M 
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Figure  13.  Waveform  L2  (0.05-Second Segment) 
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Figure  14.  Waveform  L?  (0.01-Second  Segment) 


Figure  15.  Spectrum  of  Waveform  L / 
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Figure  / 6.  Spectrum  of  Waveform  Li 
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Figure  1 7.  Spectrum  of  Waveform  M 
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7.2  FIRST-ORDER  FIT 


A  combined  least  squares  first-  and  second-order  fit  was  conducted  using  laser  L2  as 
excitation  x(n)  and  microphone  M  as  nonlinearity  response  z(n).  The  optimum  first-order 
kernel  estimate  /i,  is  displayed  in  figure  18.  The  memory  length  K  was  taken  as  336,  giving  a 
memory  duration  of  3.5  ms  to  both  the  first-  and  second-order  kernels.  This  duration  is 
apparently  long  enough  to  encompass  the  significant  contributions  of  these  kernels,  because  ht 
has  decayed  by  both  ends  of  the  fitting  interval. 

The  corresponding  estimate  of  the  first-order  frequency-domain  kernel  H]  is  shown  in 
figure  19.  The  fitting  procedure  was  limited  to  32  kHz;  the  oscillations  beyond  this  frequency 
are  irrelevant  side  lobes.  The  major  first-order  response  occurs  between  frequencies  4  kHz  and  9 
kHz.  Thus,  the  high-frequency  content  observed  in  response  z(n)  in  figure  17  is  not  due  to  a 
linear  effect  on  input  x(n),  but  is  likely  due  to  some  nonlinear,  higher  order  effect  within  the 
cicada’s  body  itself. 

7.3  SECOND-ORDER  IDEALIZATIONS 

To  more  easily  understand  some  of  the  second-order  results  to  be  presented  for  this  cicada 
investigation,  it  is  worthwhile  to  first  consider  some  idealizations  of  the  second-order  kernels, 
both  in  the  two-dimensional  time-delay  domain  r,,  r2  and  in  the  two-dimensional  frequency 
domain  /, ,/, .  The  examples  will  be  taken  in  the  continuous  time-delay  domain  for  simplicity; 
these  cases  have  their  direct  counterparts  in  the  discrete  time-delay  domain.  The  first  case  is 

=  -r,), 

(71) 

where  vv  is  a  broad  function,  and  W  is  its  Fourier  transform.  This  second-order  time-domain 
kernel  is  a  line-impulse  function  along  the  45°  line  in  r,,r2 ,  with  a  slowly  varying  envelope  w. 
The  second-order  frequency-domain  kernel  is  a  narrow  mountain  ridge  centered  along  the  -45° 
line  in  its  domain.  The  height  of  the  ridge  is  constant  along  any  -45°  line.  Both  functions  in 
equation  (71)  are  obviously  symmetric  in  their  respective  domains.  The  output  time  waveform 
corresponding  to  equation  (71)  is 


(72) 


There  is  memory  involved  in  this  squaring  device,  but  there  are  no  “cross-terms”  in  which  two 
delayed  values  ofx  are  multiplied  together,  such  as  in  equations  (1)  and  (24). 
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Figure  18.  First-Order  Kernel  hi(k),  \hj(k)\  (:),  and  Envelope  (r) 


0  I - 1 — - t - 1 - 1 - r - — i - t - 1 - r 


frequency  (kHz) 


Figure  19.  First-Order  Kernel  Ht 
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The  second  case  is 


/72(rpr2)  =  n|  -  ^  s(r,  -r2), 


H2{f,J2)  =  W(f,+f2)S 


/,-/2 


(73) 


The  function  s-  is  a  real,  even,  narrow  function,  in  which  case  its  Fourier  transform  S  is  a  real, 
even,  broad  function.  Function  w  is  arbitrary  but  generally  broad.  Time-domain  kernel  h2  has  a 
ridge  concentrated  near  the  45°  line  but  has  a  slowly  varying  amplitude  governed  by  vv. 
Frequency-domain  kernel  //,  is  concentrated  near  the  -45°  line  and  has  a  slowly  varying 
amplitude  governed  by  S.  The  time-domain  waveform  corresponding  to  kernel  (73)  is 


y2( 0  =  jj drt  Jt2 


TJ  +  T_2  ' 

2  J 


s( r,  -T2)  x(t  -Ti)  x(t-T2), 


(74) 


which  has  limited-extent  cross  terms,  governed  by  the  duration  of  s. 
The  third  case  is  given  by 


M*i>r2)  =  2  cos[2;r  fL  (r,  +t2)  +  J3]  u 


r,  +t2 


■s(*j  -t2). 


HA/, .A)  -  [exp Ofi)  W{f  +  A  -  A)  +  exp(-;«  W{A+A+A )]S 


A- A 


(75) 


Frequency  fL  causes  the  ridge  of  IT  to  move  up  and  down  in  the  fy,f2  plane  by  the  value  of  ft 
and  thereby  creates  two  ridges  parallel  to  the  ^45°  line.  If  s(0)  =  0,  then  the  time-domain  kernel 
/i2(r,,r2)  is  zero  for  r2  =  r,,  which  is  the  central  +45°  line.  Combinations  of  different 
functions  vv  and  .v,  along  with  different  parameter  values  for  f,  and  (i  lead  to  a  wide  variety  of 
possibilities  for  the  two  second-order  kernels  h2  and  H2. 

The  estimate  of  the  second-order  time-domain  kernel  h2  for  laser  L2  of  cicada  song  S14  is 
displayed  in  figures  20  and  21.  Since  this  kernel  has  a  very-high-frequency  behavior,  it 
oscillates  considerably  in  r,  and  r, .  Accordingly,  it  has  been  rectified  and  smoothed  so  that  it 
is  easier  to  observe  where  its  major  contributions  occur.  Figure  20  is  for  an  elevation 
observation  angle  of  35°,  while  figure  21  is  a  top  (90°  elevation)  view.  There  is  a  pronounced 
contribution  along  two  symmetrically  displaced  lines  at  +45°  and  a  deep  null  along  the  main 
+45°  line.  These  behaviors  are  consistent  with  those  anticipated  by  equations  (71),  (73),  and 
(75).  The  locations  of  the  main  energy  are  near  the  lines  where  |£2  -  kt  |  =  220.  At  the  sampling 

frequency  of  96  kHz,  this  corresponds  to  a  time  delay  of  2.3  ms  in  the  cross  terms  of  y2(n). 
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Figure  21.  Smoothed  Second-Order  Kernel  \lt2\  (Observation  Elevation  Angle  90°) 
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The  estimate  of  the  second-order  frequency-domain  kernel  H2  for  laser  L2  is  displayed  in 
figures  22  and  23,  the  former  for  an  observation  elevation  angle  of  50°  and  the  latter  for  an 
elevation  angle  of  90°  (top  view).  The  strong  response  (red  line  in  figure  23)  centered  at 
fi+  f2  =  0  corresponds  to  a  strong,  low-frequency  component  in  second-order  model  response 
v2(/7),  while  the  two  strong  spread  lines  centered  at  _/j  +/2  =  ±  6  kHz  correspond  to  the  spectral 
hump  for  microphone  waveform  z(n)  in  figure  17.  It  can  be  seen  from  figure  23  that  the 
second-order  frequency-domain  kernel  H2(f\  ,/2)  responds  very  well  all  the  way  up  to  32  kHz, 
which  was  the  highest  frequency  fitted.  But  that  large  response  only  takes  place  in  the  second 
and  fourth  quadrants  of  fx  ,/2  space,  which  is  the  region  where  difference  frequencies  are 
created.  There  is  very  little  response  in  the  first  and  third  quadrants;  note  that  figures  22  and  23 
are  in  dB.  (The  horizontal  and  vertical  white  streaks  in  figure  23  above  32  kHz  should  be 
ignored;  they  are  due  to  a  malfunction  in  producing  the  MATLAB  .tif  plots  for  inclusion  in  a 
Microsoft  Word  document.  Several  trials  yielded  these  same  results  with  streaks.) 

It  is  not  valid  to  compare  the  level  of  h2  with  the  level  of  /?, ;  these  two  time-domain  kernels 
don’t  even  have  the  same  dimensions;  this  can  be  seen  from  equation  (1),  where  an  extra* 
multiplication  takes  place  at  second  order.  The  same  comment  applies  to  attempting  to  compare 
the  level  of  frequency-domain  kernel  H2  with  the  level  of  Hx\  these  two  quantities  don’t  have 
the  same  dimensions  either.  The  only  valid  comparison  of  levels  is  in  terms  of  the  Volterra 
model  output  waveforms  yx{n)  and  y2(t 7),  which  always  do  have  the  same  dimensions  and 
which  are  the  same  as  those  of  z(n). 

A  plot  of  a  section  of  the  fitted  first-order  waveform  yx  (/?)  is  given  in  figure  24.  It  has  no 
jagged  behavior,  as  was  present  in  microphone  output  z(n )  in  figure  12  or  in  laser  input  x(n)  in 
figures  1 3  and  1 4.  This  is  consistent  with  the  estimated  first-order  frequency-domain  kernel  Hx 
in  figure  1 9,  which  has  its  major  response  below  1 5  kHz. 

The  corresponding  spectrum  )j(/)  of  yx(n)  is  displayed  in  figure  25.  Its  major 
contribution  is  in  the  frequency  band  from  5  to  9  kHz;  this  compares  favorably  with  the  result  for 
the  first-order  frequency-domain  kernel  //,  in  figure  19.  It  should  be  noted  in  figure  25  that  the 
zero-frequency  level  is  approximately  the  same  as  the  major  response  level  at  8  kHz.  In  contrast, 
in  figure  19,  the  zero-frequency  level  is  significantly  below  that  at  8  kHz.  This  difference  is  due 
to  the  fact  that  the  excitation  x(n)  of  laser  L2  has  a  much  higher  spectral  level  near  zero 
frequency  than  at  8  kHz,  as  can  be  seen  in  figure  16. 
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Figure  22.  Second-Order  Kerne l  H 2  (Observation  Elevation  Angle  50°) 
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Figure  23.  Second-Order  Kernel  (Observation  Elevation  Angle  90°) 
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Figure  24.  First-Order  Waveform  }’/ 


10  i - 1 - 1 - 1 - 1 - 1 - 1  i  i - r 


1 1 _ i _ I _ i _ i _ i _ i _ I _ i _ L 

0  5  10  15  20  25  30  35  40  45 


frequency  (kHz) 


Figure  25.  First-Order  Spectrum  Y/ 
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A  section  of  second-order  Volterra  component  waveform  y2(n)  is  illustrated  in  figure  26. 

It  has  a  somewhat  more  jagged  behavior  than  first-order  waveform  _y,(/ ?)  in  figure  24.  The 
corresponding  spectrum  Y2(f)  of  y2(n)  is  given  in  figure  27.  A  pronounced  spectral 
contribution  of  Y2  is  centered  at  6  kHz  with  a  peak  value  of  7  dB.  The  first-order  spectrum  7,  in 
figure  25  has  a  peak  value  of  only  -3  dB.  This  10-dB  difference  is  due  to  the  second-order 
frequency-domain  kernel  results  for  H2  in  figures  22  and  23.  Namely,  it  is  the  ±6-kHz  ridges 
of  H2  in  figure  23  that  are  creating  difference  frequencies  in  the  neighborhood  of  6  kHz.  Here, 
it  is  valid  to  make  these  quantitative  comparisons  of  levels  between  Y}  and  K,,  because  the  same 
spectral  estimation  techniques  and  parameters  were  used  on  v,(w)  and  y2(n),  which  have  the 
same  dimensions. 

The  corresponding  section  of  total  Volterra  fit  y(n)  is  plotted  (in  red)  in  figure  28, 
superposed  on  the  measured  microphone  response  z(n)  (in  black).  As  can  be  seen,  there  are 
time  segments  where  the  total  fit  is  especially  good,  and  other  time  segments  where  the  fit  is 
poor.  This  is  partially  attributed  to  the  fact  that  only  laser  L2  was  used  for  this  least  squares  fit. 
There  is  the  additional  input  from  laser  L,  that  was  totally  ignored  during  this  single-input  least 
squares  fit. 

All  the  fitting  results  thus  far  have  been  for  laser  L,  of  cicada  song  S14  considered  as  the 
input.  For  the  fitting  results  in  figure  29,  laser  L,  was  used  instead  as  the  input  x(n).  The 
measured  microphone  output  z(n)  is  plotted  in  black,  fit  y{(n)  is  plotted  in  blue,  and  fit  y2(n) 
is  plotted  in  red.  For  this  segment  of  time,  the  first-order  fit  is  extremely  small,  while  the 
second-order  fit  does  a  decent  job  of  mimicking  output  z(n).  In  fact,  for  this  particular  time 
segment,  the  ratio  of  powers  in  >’,(»)  to  z(n)  is  0.036,  while  the  ratio  of  powers  in  y2(n)  to 
z(n)  is  0.68.  This  example,  along  with  figure  28,  brings  out  the  point  that  the  cicada  with  two 
laser  measurements  and  one  microphone  measurement  should  be  considered  as  a  /vvo-input,  onc- 
output  nonlinear  system  with  memory.  That  is,  the  block  diagram  in  figure  1  is  not  general 
enough  to  encompass  this  generalization.  More  will  be  said  on  this  observation  below. 

Results  for  a  different  cicada  example  (namely,  song  S9)  are  given  in  figures  30  and  3 1  for 
lasers  L,  and  L-, ,  respectively.  Again,  the  measured  microphone  output  z(n)  is  plotted  in 
black,  fit  v,  (/?)  is  plotted  in  blue,  and  fit  y2(n)  is  plotted  in  red.  The  fractional  power  fit  of  total 
waveform  y(n)  to  output  z{n)  in  figure  30  for  laser  L,  is  0.25  at  the  left  edge  of  the  figure,  but 
0.92  at  the  right  edge  of  the  figure.  By  contrast,  the  fractional  power  fit  of  total  waveform  y(n) 
to  output  z(n)  in  figure  3 1  for  laser  L2  is  0.78  at  the  left  edge  of  the  figure,  but  0.24  at  the  right 
edge  of  the  figure.  This  example  serves  to  further  support  the  fact  that  different  tymbals  (lasers) 
are  contributing  to  different  portions  of  the  time  waveform  emitted  from  the  cicada.  The  only 
way  to  properly  calculate  the  relative  strengths  of  these  components  is  by  way  of  a  two-input, 
one-output  Volterra  model. 


47 


x  10‘J 


Figure  26.  Second-Order  Waveform  y 2 


Figure  27.  Second-Order  Spectrum  Y2 
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Figure  28.  Total  Waveforms  z(k)  y(r) 


x  10'J 


Figure  29.  Waveforms  z(k )  yfb)  y2(r )  (Laser  L2  of  Cicada  Song  SI  4) 
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Figure  30.  Waveforms  z(k)  yi(b)  )'y(r)  (Laser  L /  of  Cicada  Song  S9) 
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Figure  31.  Waveforms  z(k)  yt(b)  yi(r)  (Laser  L2  of  Cicada  Song  S9) 
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8.  SECOND-ORDER,  TWO-INPUT,  ONE-OUTPUT  VOLTERRA  MODEL 


A  Volterra  model  of  a  second-order,  two-input,  one-output  system  is  given  in  figure  32.  The 
two  inputs  are  now  labeled  x,  ( n )  and  x2(n).  There  are  now  five  required  kernels  that  must  all 
be  estimated  from  simultaneous  recordings  of  laser  inputs  x,(«)  and  x2(/?)  and  microphone 
response  z(n).  The  number  of  subscripts  on  the  kernels  h  denotes  their  order,  either  first  order 
or  second  order.  A  new  type  of  kernel  must  now  be  employed,  namely,  a  “cross-kernel'’ 
h]2(k\,k2 ) ,  which  has  two  inputs  and  one  output  and  which  yields  a  second-order  waveform 


A.  I  A'  I 

.V12 (^) =  T.  ^12(^1  *^2 )  •*)  )  x2 (n  —  k2 ) 

k  =0  *,=0 


=  A:  \\df,  df2  cxp[/2^(/,  +/2)/;A] //,,(/,,/,)  A, (/,)  X2(f2). 


(76) 


However,  time-domain  cross-kernel  hl2(kt,k2 )  need  not  be  symmetric,  because  two  different 
inputs  x  are  multiplied  in  the  first  line  of  equation  (76).  Therefore,  the  only  symmetry  relation 
satisfied  by  its  corresponding  second-order  frequency-domain  kernel  is 


(77) 


which  follows  from  the  imposed  realness  of  /?,,(£,,£,). 

If  Volterra  output  yl2(n)  is  to  have  frequency  content  only  in  the  band  for 

purposes  of  fitting  to  a  band-limited  version  of  z(n),  and  if  Xt(f)  and  X2(f)  are  broadband, 
the  second  line  of  equation  (76)  reveals  that  one  must  have  frequency-domain  kernel  H]2  (/,  ,f2 ) 
nonzero  only  for 

<|./i+/2|</^  m 

just  as  in  equation  (34).  Therefore,  it  will  be  necessary  to  let  variables  /,,/,  in  figure  2  range 
not  only  over  the  blue  region,  but,  in  addition,  over  its  green  extension  into  the  second  quadrant. 
The  number  of  unknown  (complex)  coefficients  for  //,,  w  ill  be  equal  to  the  sum  of  the  number 
of  coefficients  required  for  each  of  //,,  and  H 22 .  This  larger  number  of  coefficients  will 
significantly  increase  the  effort  required  in  the  least  squares  procedure. 

The  basis  functions  for  the  /7I2(A',A,A',A)  kernel  are  given  by  one-half  of  equation  (36)  and 
its  conjugate  (but  not  using  the  symmetry  property  in  the  top  line  of  equation  (37)),  namely, 

a{mx,m2)  cos(/l)  +  b(mt,m2 )  sin(/l),  (79) 
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where  A  is  given  by  equation  (38).  The  b(mx  -mx )  terms  must  again  be  dropped  from  the  basis 
set. 


h 


Figure  32.  Second-Order,  Two-Input,  One-Output  Volterra  Model 
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9.  SUMMARY 


Several  new  concepts  have  been  introduced  to  achieve  some  alleviation  of  the  curse  of 
dimensionality  in  determining  the  higher  order  kernels  in  a  Volterra  expansion.  These  include: 
partitioning  of  the  frequency  scale(s),  the  fundamental  region  of  formation,  the  diagonal  strips  in 
the  two-dimensional  frequency  space,  the  second-  and  third-order  basis  functions,  and  filtering  of 
the  measured  nonlinear  output  to  the  current  frequency  band  under  investigation.  The  equations 
for  the  real  basis  functions  at  second  and  third  order  take  some  unexpected  forms  that  could 
probably  not  be  anticipated.  These  results  were  based  on  taking  full  advantage  of  symmetry  and 
conjugate-symmetry  frequency-domain  relations  that  exist  for  real,  symmetric,  second-order  and 
third-order,  time-domain  kernels.  The  inability  to  construct  ideal  bandpass  filters  requires  the 
use  of  a  frequency-overlap  procedure,  follow  ed  by  the  discarding  of  edge  estimates  with  inherent 
errors,  and  the  retention  of  only  the  interior  estimates  of  higher  accuracy. 

Application  to  a  first-  and  second-order  (noise-free)  control  example  with  a  white  broadband 
excitation  gave  excellent  estimates  of  all  the  first-  and  second-order  properties,  such  as  the 
individual  kernels,  the  individual  Volterra  waveforms,  and  the  total  estimated  output  waveform. 
Application  to  a  cicada  mating  call  with  a  distinctly  non-white  and  non-Gaussian  excitation  gave 
good  results  for  the  estimated  first-  and  second-order  kernels  and  waveforms,  considering  the 
non-optimality  of  this  type  of  excitation.  However,  if  the  design  of  the  input  excitation  to  a 
nonlinear  system  with  memory  is  under  the  user’s  control,  the  ideal  excitation  to  use  is  a  white, 
Gaussian  process.  Better  correlation  properties  ensue  and  smaller  condition  numbers  prevail  in 
this  latter  situation. 

The  reduction  factor  in  the  number  of  coefficients  to  be  determined  at  each  frequency 
partition  is  given  by  W/F,  w  here  W  is  the  bandwidth  of  the  frequency  interv  al  under 
investigation,  and  F  is  the  largest  frequency  of  interest,  typically  near  the  Nyquist  frequency. 

This  factor  applies  at  all  orders  and  is  helpful  in  alleviating  the  curse  of  dimensionality,  but  does 
not  eliminate  it. 
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